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Chapter  1 

Executive  Summary 


1.1  Thrust 

Phillips  Laboratory  has  recently  realized  the  confluence  of  two  important  events:  the  dis¬ 
covery  of  a  spatially  discreet  transport  theory  and  the  construction  of  a  cellular  automata 
machine  (CAM-8)  and  a  wide-purpose  Connection  Machine-5  (CM-5).  That  is,  there  now 
exists  first-principles  lattice  gas  automata  (LGA)  and  lattice  Boltzmann  equation  (LBE) 
formalisms  for  modeling  complex  systems.  The  CAM-8  now  exists  as  a  cheap,  fast  paral¬ 
lel  bit-level  hardware  optimized  for  LGA  simulation.  The  Connection  Machine-5,  which  has 
been  recently  installed  at  the  Army  High  Performance  Computing  Research  Center,  is  ideally 
suited  to  LBE  simulation  because  of  its  floating-point  and  virtualization  capabilities.  There¬ 
fore,  we  propose  a  two-pronged  parallel  computing  strategy  for  our  thermohydrodynamic 
research:  1)  LGA  implemented  on  a  low-cost  next-generation  cellular  automata  machine 
(CAM-8);  and  2)  LBE  implemented  on  the  Connection  Machine  (CM-200  and  CM-5). 

1.2  Approach 

1.2.1  Lattice  Gas  Automata 

LGA  and  LBE  offer  a  simplified  modeling  strategy  for  handling  complex  and  “messy”  physi¬ 
cal  simulations,  for  example,  fluids  with  complex  interfacial  boundaries  throughout  the  entire 
simulation  space.  LGA  and  LBE  are  an  efficient  software  tool  for  programming  massively 
parallel  architectures.  LGA  offers  a  unique  and  powerful  representation  of  macroscopic  dy¬ 
namics  by  reducing  the  calculation  to  local  space-time  processes  controlled  within  standard 
dynamic  random  access  memory  employing  only  data  shifts  and  permutations.  No  central 
processing  unit  operations  are  required  The  payoff  of  a  widely  successful  LGA  theory,  can 

^The  CAM-8  design  requires  no  internal  CPUs. 
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be  high. 

CAM  is  a  promising  platform,  comparable  in  speed  ^  and  less  expensive  then  other 
parallels  platforms  in  the  United  States  capable  of  running  LGA  calculations  We  will  test 
the  validly  of  massively  hne-grained  simulation  and  its  correspondence  with  hydrodynamic 
macroscopic  systems. 

1.2.2  Vision 

We  wish  to  exploit  LGA  and  LBE  theory  and  next  generation  massively  parallel  architectures 
for  boasting  simulation  rates  and  cost-effectiveness.  Potential  long-range  application  is: 
accurately  model  the  time  evolution  of  complex  Navier-Stokes  fluids. 

1.3  Anticipated  Benefits 

1.3.1  Need  of  the  Simulation  Community 

This  initiative  element  supports  the  Air  Force  critical  technology  of  simulation  and  modeling, 
an  important  new  growth  area  for  our  laboratories.  The  Air  Force  objective  for  broad-based 
environmental  simulations  to  help  reduce  future  design  and  engineering  costs  are  critically 
dependent  on  physical  simulation  technology.  A  sufficiently  good  physical  simulation  strategy 
and  platform  does  not  yet  exist.  The  direction  in  which  the  Air  Force  simulation  community 
is  heading,  although  highly  appealing  to  Air  Force  leaders,  is  actually  technically  very  risky. 
Basic  research  in  new  parallel  computing  strategies  is  needed  to  help  reduce  that  risk. 


^Data  shifting  is  handled  simply  by  pointer  manipulations. 

^MasPar’s  MP-1  is  a  typical  example  with  8192  four-bit  processors  costing  about  $800K. 


Chapter  2 

Programmatic  Description 

2.1  Introduction 

Request  AFOSR  support  for  a  new  initiative  for  lattice  gas  automata  research  and  develop¬ 
ment  of  a  massively  parallel  cellular  automata  machine  (CAM).  This  research  will  comprise 
constructing  a  large  CAM  and  developing  hydrodynamic  lattice  gas  automata  (LGA)  algo¬ 
rithms.  This  research  will  also  comprise  implementing  the  lattice  Boltzmann  equation  (LBE) 
on  state-of-the-art  parallel  supercomputer  architectures  such  as  the  Connection  Machine 
(CM-200  and  CM-5).  The  hydrodynamic  LGA  will  serve  as  a  test-bed  for  demonstrating 
a  novel  parallel  computing  strategy  which  strives  to  capture,  in  the  most  digitally  efficient 
and  numerically  stable  way,  accurate  physical  behavior  of  large-scale  complex  systems.  This 
research  will  also  provide  for  comparison  of  the  LGA  method  with  LBE.  We  will  explore 
using  a  multigrid  approach  in  connection  with  the  lattice  Boltzmann  method. 

2.2  Submitting  Agency 

Phillips  Laboratory  (Air  Force  Materiel  Command),  Geophysics  Directorate,  (PL/GPAA, 
Hanscom  AFB,  MA  01731).  Manager;  J.  Yepez,  DSN  478-2475,  yepez@plh.af.mil.  TAP 
Reference:  Gomputational  Mathematics  Subarea  2304/GS,  PL/GP  Simulation  and  Applica¬ 
tions  6670. 

2.3  Lattice-Gas  Automata  and  Lattice  Boltzmann  Model 

LGA  and  LBE  methods  are  still  in  a  formative  stage.  Important  engineering  applications 
for  hydrodynamic  LBE  currently  only  exist  for  flow  through  porous  media.  LBE  is  suited 
to  parallel  architectures  with  significant  floating-point  performance,  high  virtual-processor 
ratios,  and  efficient  grid-communications.  The  Gonnection  Machine  is  an  ideal  platform 
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for  LBE.  On  the  other  hand,  LGA  is  suited  to  parallel  architectures  optimized  for  hne- 
grained  bit-level  manipulations.  Cellular  automata  machines  (CAM)  can  be  constructed 
inexpensively  to  run  LCA  models  at  unprecented  rates,  orders  of  magnitude  faster  than 
LCA  models  on  general-purpose  massively  parallel  computers. 

A  3-D  hydrodynamics  model  is  a  prototypical  complex  system  for  the  LCA  or  LBE 
methodology.  LCA  and  LBE  are  an  efficient  and  physically  elegant  parallel  programming 
paradigm.  All  calculations  are  local  and,  for  LCA,  can  use  reversible  rules.  Therefore 
the  LCA  and  LBE  programming  paradigms  closely  mirror  real-world  physics  which  is  hue- 
grained  and  time-invariant.  So  LCA  and  LBE  may  offer  a  simplihed  modeling  strategy 
for  handling  many  complex  and  “messy”  physical  simulations,  i.e.  cases  with  multiphase 
interfacial  boundaries  affect  huid  how  properties.  Hydrodynamics  LCA  and  LBE  may  prove 
to  be  the  most  efficient  software  tool  for  programming  massively  parallel  architectures. 

2.4  Value  to  the  Scientific  Community 

Cellular  automata  provide  a  conceptually  simple  and  intuitive  model  for  encoding  the  essen¬ 
tial  physics  of  complex  systems.  Cellular  automata  models  of  such  systems  oher  an  unique 
opportunity  to  explore  a  system’s  underlying  statistical  behavior.  Probability  distributions, 
correlation  functions,  transport  properties,  etc.,  can  be  explored  in  detail.  LCA,  as  a  discreet 
form  of  molecular  dynamics  simulation,  provides  the  researcher  with  complete  system  de¬ 
tails  not  obtainable  by  either  empirical  or  analytical  treatments.  LCA  as  a  generalization  of 
cellular  automata,  ohers  a  unique  and  powerful  representation  of  macroscopic  dynamics  by 
reducing  the  calculation  to  local  space-time  processes,  collisions  and  streaming,  controlled 
within  standard  dynamic  random  access  memory  by  employing  only  shifts  and  permuta¬ 
tions  of  cell  data  LCA  models  conserve  exactly,  throughout  the  entire  simulation  run, 
all  physical  moments  of  the  system,  i.e.  mass,  energy,  linear  momentum,  angular  momen¬ 
tum.  Consequently,  LCA  may  prove  to  be  a  numerically  stable  modeling  technique  with 
performance  and  accuracy  approaching  hnite-diherence  methods. 

LBE  is  a  generalization  of  LCA,  where  single-particle  distributions  are  encoded  directly 
using  real  numbers.  The  simulation  dynamics  is  driven  toward  an  equilibrium  distribution 
which  the  modeler  analytically  determines  a  prior  and  encodes  into  the  simulation  algorithm. 
The  beneht  of  LBE  over  LCA  is  noise  reduction.  However,  this  beneht  is  gained  at  the  cost 
of  giving  up  exact  conservation  of  mass,  energy,  an  momentum  for  a  statistical  conservation 
of  these  quantities  and  also  at  the  cost  of  requiring  expensive  float-point  hardware.  We 

^Update  of  cell  data  can  be  achieved  by  using  fast  look-up  tables.  In  principle,  no  central  processing  unit 
operations  are  required. 
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understand  that  the  actual  benefit  of  LGA  and  LBE  methods  over  more  traditional  methods 
is  clearly  uncertain.  However,  LGA  and  LBE  merit  the  investment  of  basic  research  time 
and  resources  because  the  potential  payoff  of  such  methods,  if  successful,  can  be  very  large. 


2.5  Anticipated  Benefits  to  the  US  Air  Force 

Recently  simulation  technology  has  gained  considerable  recognition  at  the  highest  levels  in 
the  Air  Force.  Simulation  technology  is  now  considered  a  critical  technology  and  an  im¬ 
portant  new  growth  area  for  our  laboratories.  At  Phillips  Laboratory  we  have  established 
the  Office  of  Environmental  Simulation.  This  office  will  help  coalesce  DoD  resources  re¬ 
lated  to  environmental  simulation  into  a  state-of-the-art  technology  base  from  which  system 
program  offices,  DoD  engineering  and  research  components,  major  defense  contractors,  etc., 
in  the  future  may  draw  upon  as  they  design  new  Air  Force  systems  which  are  impacted 
by  the  environmentAn  particular,  tropospheric  effects.  The  capability  to  do  wide  ranging 
computer-based  proof-of-concept  tests  of  new  system  designs  using  simulation  technology  is 
hoped  to  be  a  cost  effective  measure  leading  toward  reduced  acquisition  costs.  This  is  per¬ 
ceived  by  Air  Force  leaders  as  an  critical  technology  shift  to  stave  off  adverse  impacts  of  the 
declining  defense  department  bndgets.  Yet  this  direction  in  which  the  Air  Force  community 
is  heading,  althongh  highly  appealing  to  most  Air  Force  leaders,  is  actually  technically  very 
risky. 

Althongh  the  objective  is  for  broad-based  environmental  simulations  to  reduce  future 
design  and  engineering  costs,  this  has  in  no  way  been  proven  a  viable  ronte.  Gnrrently  6.2 
based  simulation  technology  thrnsts  are  driven  by  non-physical  constraints.  Often  a  model’s 
visualization  and  graphics  output  is  viewed  as  the  critical  ingredient  to  “good”  simulation. 
Graphics  and  visnalization  are  easily  nnderstood  and  conseqnently  seem  to  receive  high  pri¬ 
ority.  Underlying  physical  components  of  snch  models  are  either  missing  or  underemphasized 
and,  therefore,  in  serious  risk  of  large  inaccuracy  and  numerically  unstable  behavior.  Basic 
research  using  efficient  or  “qnick”  methods  to  explore  certain  physical  properties  of  an  en¬ 
vironmental  system  will  help  to  underprop  the  Air  Force  simulation  thrust.  LGA  is  a  novel 
method  for  doing  such  quick  physical  explorations. 

Most  simulation  technology  project  managers  would  like  robust  and  physically  based 
engines.  But  these,  for  the  most  part,  do  not  exist  because  on  the  one  hand  there  is  the 
nearly  overwhelming  complexity  of  the  problems  and  on  the  other  hand  there  is  the  limited 
compntational  performance  of  desktop  workstation  class  machines  on  which  most  current 
day  simulation  models  reside. 

The  ambitious  goals  of  the  simulation  commnnity  are  not  synchronized  with  the  latest 
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supercomputer  technology.  This  must  be  corrected  for  any  real  progress  to  be  made  in  this 
held.  The  principle  platform  for  current  Air  Force  simulation  technology  is  the  workstation. 
Although  the  workstation’s  performance  has  reached  the  point  where  traditional  mainframes 
are  perhaps  no  longer  as  practical  as  they  once  were,  workstation  performance  is  still  orders 
of  magnitude  below  supercomputing  levels. 

The  principle  objective  of  this  basic  research  initiative  is  to  explore  a  novel  parallel  com¬ 
puting  strategy  which,  if  successful,  will  serve  as  an  additional  tool  to  aid  the  Air  Force  sim¬ 
ulation  thrust.  We  propose  building  a  large  hne-grained  cellular  automata  machine  (CAM), 
still  small  enough  for  office  desk-side  use,  yet  with  a  performance  comparable  to  a  mas¬ 
sively  parallel  Connection  Machine-200.  I  do  not  believe  that  such  a  CAM  will  solve  the 
simulation  grand  strategy  envisioned  by  Air  Force  leaders,  certainly  not  in  the  near  future. 
Nor  do  I  believe  that  such  a  CAM  will  serve  as  general-purpose  computer  for  engineering 
tasks  and  scientihc  exploration,  again,  certainly  not  in  the  near  future.  But  I  do  believe 
that  such  a  CAM  is  a  promising  platform,  and  being  optimized  for  lattice  gas  automata 
(LGA)  calculations,  is  an  ideal  platform  for  doing  quick  physical  explorations  of  complex 
fluid  systems. 

The  LGA  approach,  which  can  be  viewed  as  a  discreet  molecular  dynamics  (MD)  model¬ 
ing  approach  over  a  discreet  space  and  time,  and  the  LBE  approach,  which  can  be  viewed  as 
a  statistical  discreet  MD,  are  known  to  have  application  to  hydrodynamics,  thermohydrody¬ 
namics,  magnetohydrodynamics,  reaction-diffusion  systems,  polymer  melts.  In  the  future, 
it  is  likely,  the  LGA  and  LBE  approaches  should  have  application  to  geophysical  problems 
involving  hne-grained  locally  interacting  dynamics.  With  this  new  initiative  we  propose 
building  a  large  cellular  automata  machine  and  doing  a  detailed  study  of  three-dimensional 
hydrodynamics  including  multiphase  huids.  Our  purpose  is  to  see  if  large  CAM  machines 
can  be  constructed  and  can  handle  complex  huid  modeling  problems. 

Although  our  hydrodynamics  work  may  prove  to  be  an  effective  tool  in  certain  messy 
regimes  where  traditional  methods  are  not  well  suited,  this  would  be  only  one  of  the  pay-offs 
of  our  new  initiative  research.  Another  pay-off  will  be  a  proof-of-concept  of  a  novel  parallel 
computing  strategy  which  captures  in  the  most  digitally  efficient  and  numerically  stable  way, 
the  accurate  physical  behavior  of  complex  huid  systems. 

2.6  Funding 

The  projected  total  cost  for  the  proposed  basic  research  is  supported  by  the  Air  Force  Office 
of  Scientihc  Research  contributing  core  funding  and  new  initiative  funding  over  the  hve-year 
project  cycle,  hscal  year  1994  to  1998.  Fiscal  year  1994  funding  is  seed  funding.  The  Phillips 
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Laboratory  expenses  include  salaries  for  in-house  research  scientists  and  engineers  as  well 
as  CAM  design  and  construction  costs.  Out  year  hgures  for  projected  6.2  and  6.3  funding 
level  transition  have  not  been  estimated.  Also  not  estimated  are  supercomputer  service 
unit  costs  for  the  Connection  Machine  200  and  Connection  Machine  5  at  the  Army  High 
Performance  Computing  Research  Center  (AHPCRC).  Resources  at  AHPCRC  are  requested 
under  a  separate  but  complementary  proposal. 

2.7  Tasks 

•  Develop  3-D  hydrodynamic  lattice  gas  and  lattice  Boltzmann  automata  code  including 
multiphase  fluid  system.  Explore  lattice  Boltzmann  dynamics  which  can  model  a  fluid 
system  with  a  general  equation  of  state. 

•  Test  novel  lattice  gas  computing  architectures:  128  million  16-bit  site  CAM-8  pro¬ 
totype.  Design  1-billion  16-bit  site  CAM-8-64  machine  64  modules  using  an  SBIR 
vehicle. 

•  Explore  ways  of  using  multgrid  methods  in  lattice  Boltzmann  simulations. 

2.8  Schedule  and  Milestones 


5 


FY94 


FY95 


FY96 


TASKl 

Develop  3D  LGA  on  CAM-8 
and  LBE  on  CM-5 

LBE  with  General  Equation  of  State 

Statistical  Studies  of  Transport  Properties 

Multiphase  Fluid  with  Latent  Heating 

Dynamics  Studies:  3D  Convection, 
Thermal  Topography,  Radiation-Transfer 

TASK  2 

128-Module  CAM-8  Design 
Design  Review 
Assembly  and  Testing 
Operational  Unit 

TASKS 

Multigrid  Methods  for  LBE 


Figure  2.1:  Schedule  and  Milestones 
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Chapter  3 


Hydrodynamic  Lattice-Gas 
Automata 

3.1  Background 

In  1983  Stephan  Wolfram  proposed  a  simple  mathematical  model  to  investigate  self-organization 
in  statistical  mechanics[3].  He  proposed  a  new  “elementary”  formalism  which  he  termed  cel¬ 
lular  automata  (CA)  to  treat  dynamical  systems  on  a  discreet  spatial  and  temporal  lattice 
where  each  site  or  “node”  of  the  lattice  has  a  binary  0  or  1  value.  The  conhguration  of  the 
discreet  lattice  then  records  the  state  of  the  dynamical  held.  The  evolution  of  the  held  is  de- 
hned  by  a  particular  set  of  boolean  rules.  Computational  efficiency  follows  from  this  discrete 
formalism  when  the  CA  models  a  complex  natural  system  with  a  large  number  of  simple 
identical  components  requiring  only  local  interaction.  In  this  way,  the  CA  model  may  be 
implemented  on  special  massively  parallel  computers  allowing  one  to  update  the  state  of  the 
dynamical  held  at  each  lattice  site  simultaneously.  This  becomes  computationally  efficient 
when  each  site  update  requires  knowledge  of  only  its  nearest  neighbors  and  thereby  minimizes 
grid  communications  between  processors  of  the  parallel  computer.  Wolfram  has  shown  that 
such  simple  CA’s  possess  universal  features  of  complex  nonlinear  dynamical  systems.  Since 
Wolfram’s  seminal  research,  many  applications  of  CA  have  been  found  in  diverse  helds  in¬ 
cluding  biology,  chemistry,  mathematics,  and  physics.  In  physics,  CA’s  have  had  application 
in  such  areas  as  dihusion  [9],  hydrodynamics,  information  theory,  magnetohydrodynamics 
[11],  and  magnetic  systems  [10]. 

Frisch  et.  al.  have  presented  a  CA  implementation  of  the  Navier-Stokes  equation  using  a 
lattice  gas  automata  in  a  two-dimensional  triangular  lattice  [12] .  In  two-dimensions,  for  single 
speed  automata,  only  on  a  triangular  lattice  does  the  momentum  flux  tensor  take  the  correct 
form;  and  therefore,  isotropy  is  upheld.  This  hrst  lattice  gas  automata  (LGA)  simulation 
of  two-dimensional  hydrodynamic  flow  is  call  the  FHP  model  named  after  its  authors  Frish, 
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Hasslacher,  Pomeau.  Extension  of  the  FHP  model  to  three-dimensions  proved  difficult  in 
that  momentum  conservation  was  violated,  and  isotropy  problems  again  arose.  This  was 
remedied  by  Dhumieres  et.  ah  by  employing  a  four-dimensional  face-centered  cubic  lattice 
projected  on  to  three-dimensions  to  implement  a  LGA[15].  Much  numerical  simulation  of 
2D  EGA  has  been  performed  [18]  and  comparisons  made  between  simulation  and  theory 
[19].  Therefore,  there  are  well-dehned  prescriptions  for  determining  fluid  parameters  such  as 
Rayleigh-Taylor  instability  growth  rates  [20]  and  shear  viscosity  [21,  22,  23]. 

3.2  Massively  Parallel  Implementation 

Gellular  automata  potentially  offer  a  large  beneht  over  conventional  numerical  techniques 
by  their  inherent  computational  stability  and  hne  granularity.  Yet  GA’s  promise  of  be¬ 
ing  more  “economical”  than  standard  finite-difference  schemes  used  to  numerically  solve 
partial  differential  equations  depends  upon  whether  or  not  they  can  be  rigorously  shown 
to  emulate  Navier-Stokes  dynamics.  Much  of  the  work  of  EGA  modeling  has  been  done 
on  two-dimensional  lattices  and  more  recently  on  four-dimensional  lattices  projected  on  to 
three-dimensions.  The  range  of  validity  of  the  EGA  implementations  is  still  an  important 
issue  requiring  further  evaluation.  Fundamentally,  if  we  can  validate  EGA  rules  encoding 
Navier-Stokes  dynamics  we  should  expect  a  drastic  decline  in  elapsed  computation  time  as 
we  can  spread  the  EGA  calculations  over  a  larger  and  larger  number  of  processors. 

In  concert  with  the  new  GA  techniques,  new  massively  parallel  architectures  also  abound. 
The  Gonnection  Machine  is  one  such  platform  [37]  [38].  Our  research  in  this  area  therefore 
hinges  upon  placing  a  dedicated  massively  parallel  machine  on  site.  Phillips  Eaboratory  is 
currently  evaluating  candidate  architectures  and  is  conducting  an  informal  held  survey.  The 
GM-2  and  GM-5,  currently  at  the  ARMY  High  Performance  Gomputing  Research  Genter 
(AHPGRG)  and  at  several  universities  in  Massachusetts,  is  a  convenient  platform  since  the 
G-star  and  FORTRAN  90  parallel  compilers  are  well  advanced  and  ease  the  parallel  coding 
workload.  Machines  with  the  Mach  operating  system,  which  is  based  on  message  passing, 
are  also  becoming  available.  However,  despite  the  numerous  commercial  parallel  computers 
available,  the  optimum  platform  for  EGA  simulations  are  dedicated  architectures  such  as  the 
GAM-8  currently  development  by  Margolus  and  Toffoli  at  the  MIT  Eaboratory  for  Gomputer 
Science  [40]. 


Appendix  A 

Lattice  Boltzmann  Gas  with  a 
General  Equation  of  State 

A.l  Abstract 

We  present  a  simple  way  to  add  an  arbitrary  equation  of  state  to  a  automaton  gas  modelled 
in  the  lattice  Boltzmann  limit.  As  a  way  of  interpreting  the  lattice  Boltzmann  equation  we 
present  a  new  derivation  based  on  an  automaton  Hamiltonian  and  the  Liouville  equation. 
A  convective-gradient  term  added  to  the  LBE  is  shown  to  be  a  sufficient  route  for  modeling 
hydrodynamic  flow  with  a  general  equation  of  state.  The  method  generalizes  to  multi-speed 
gases  in  three-dimensions. 

A. 2  Introduction 

Lattice  gas  methods  for  hydrodynamic  flow  over  a  discreet  hue-grained  space-time  imple¬ 
mented  on  massively  parallel  machines  like  the  CM-2  and  CM-5  [37]  or  programmable  matter 
machines  like  the  CAM-8  [40]  represent  an  important  new  avenue  for  practical  simulations  of 
complex  physical  systems.  Local  streaming  and  collision  rules  dehne  a  mesophysical  world 
underlying  the  macroscopic  system  of  interest.  This  robust  computational  methodology 
provides  an  exciting  alternative  to,  and  not  simply  an  approximation  of,  the  usual  par¬ 
tial  differential  equation  method  of  description  and  the  associated  hnite  difference  schemes. 
The  cellular  automaton  formalism,  popularized  in  the  physics  community  by  Wolfram  [4], 
has  been  extended  to  a  more  general  lattice  gas  formalism  [12,  15].  In  lattice  gas  codes, 
all  individual  boolean  bits,  representing  automaton  particles,  simultaneously  propagate  and 
rearrange  within  parallel  architectures  built  from  low-cost  digital  circuits. 

Recently,  this  innovative  lattice  gas  approach  has  been  extended  to  the  lattice  Boltzmann 
approach  [31].  In  place  of  the  detailed  streaming  and  colliding  of  digital  bits,  one  focuses 
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on  the  statistical  regime  where  only  total  particle  count  per  lattice  direction  per  lattice 
node  influences  the  dynamics,  and  so  one  neglects  particle-particle  correlations.  Although 
this  new  approach  is  less  noisy,  it  relies  on  expensive  floating  point  calculations.  The  most 
practical  simulation  method  for  production  work  may  lie  somewhere  between  the  lattice  gas 
and  lattice  Boltzmann  extremes.  Yet  the  lattice  Boltzmann  approach  offers  an  important 
analytical  advantage.  One  may  capture  the  essential  physics  of  the  complex  system  by 
stating  no  more  than  the  system’s  equilibrium  distribution. 

Here  we  exploit  the  analytical  facility  of  the  lattice  Boltzmann  approach  and  show  that 
the  addition  of  a  convective-gradient  term  in  the  lattice  Boltzmann  equation  (LBE)  allows 
one  to  model  a  hydrodynamic  gaseous  flow  governed  by  a  general  equation  of  state.  We 
restrict  ourselves  to  single-speed  automata.  Thus,  the  system  pressure  may  depend  on  local 
density  variations.  It  is  straightforward  however  to  generalize  our  result  to  a  multi-speed 
automaton  gas  so  that  the  pressure  dependence  includes  local  temperature  variations  as  well. 


A. 3  The  Lattice  Boltzmann  Equation 


We  wish  to  consider  a  simple  two-dimensional  lattice  Boltzmann  gas  defined  on  a  discreet 
spatial  lattice.  Automaton  particles  stream  through  the  lattice  and  undergo  collisions  in 
a  similar  fashion  to  the  usual  FHP  model  [12].  Therefore,  there  exists  a  small  number  of 
momentum  states  given  by 

(A.l) 


27ra  .  27ra 


ea  =  ^^cos— ,sm  — 

where  a  =  1,2, ...  ,b.  The  automaton  Hamiltonian  is  the  difference  of  kinetic  and  collision 
terms 

H  =  ^ea'  -  F  ■  q.  (A.2) 

Hamilton’s  equations  for  an  automaton  are 


and 


.  dH 

ap 


dH  ^ 

p  =  =  F. 


(9q 


(A.3) 


(A.4) 


Liouville’s  equation  for  the  distribution  function  may  be  written  in  terms  of  the  automaton 
Hamiltonian  using  the  Poisson  bracket  notation  [32]  as 


where 


\f  m  =  —  ^ 

^  5q  ■  9p  dq'  dp' 


(A.5) 

(A.6) 
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To  first  order  and  near  equilibrium  we  may  approximate  a  change  in  the  distribution  function 
by  a  collision  as 

dfa  _  (A  7\ 

dp  ||F||r  > 

where  we  have  taken  6  fa  =  fa—fa'^  and  introduced  the  collision  relaxation  time,  r,  by  writing 
=  ||-F||t.  Using  Hamilton’s  automaton  equations  (A. 3, A. 4)  and  our  linear  approximation 
of  (9/a/(9p  (A. 7),  the  Poisson  bracket  becomes 

[/a,  H]=e^-  —  + - - - ,  (A.8) 

and  so  we  obtain,  to  hrst  order,  lattice  Boltzmann  equation 

df  1 

^  +  =  (A.9) 

Shiyi  Chen  et.  ah  have  arrived  at  (A.9)  by  expanding  the  lattice  Boltzmann  collision  term 
to  hrst  order  [27]  about  an  equilibrium  distribution 

df 

^  +  ea-V/a  =  Ha.  (A.IO) 

^  =  nair)  +  7^,  (A.ll) 

where  the  zeroth-order  term,  f2a(/®'^),  must  vanish  by  construction.  The  simplest  ansatz  for 
dfla{f^‘^)/dfb  is  to  choose  it  to  be  diagonal 


dnair'i)  _  1 

dft  T 


(A.12) 


where  again  r  is  the  characteristic  relaxation  time  for  the  simulation.  Integrating  (A.12) 
leads  to  the  same  collision  term  as  in  the  R.H.S.  of  (A.9) 

n.  =  -l(/a (a.13) 

T 


The  collision  term  is  proportional  to  the  difference  of  the  distribution  function  and  its  equi¬ 
librium  value. 


A. 4  Equilibrium  Distribution 

To  compute  the  system  dynamics  on  a  parallel  machine  we  implement  the  exact  cellular  form 
of  (A.9) 

/a(x-h  ea,t-h  1)  =  /a(x,t)  -  -  (/a(x,t)  -  /a''''(x,t))  .  (A.14) 

T 


11 


The  local  cellular  automaton  rule  (A.  14)  does  not  explicitly  show  any  mixing  of  particle  flow 
directions.  That  fact  that  (A.  14)  does  represent  collisional  mixing  is  implicitly  built  into  the 
form  of  which  depends  on  the  local  density  and  flow  velocity 


P  =  H/a 

a 


(A.15) 


and 


u  = 


Ha  fa 


(A.16) 


Chen  et.  ah  [31]  have  introduced  a  pressure-corrected  lattice  Boltzmann  equation  (PCLBE) 
by  taking  the  equilibrium  distribution  to  have  the  following  Chapman-Enskog  expansion 


h‘^  =  d+^e. 


u  -h  p 


D{D  +  2) 
2c% 


('aif'ajUi'Uj 


pD  2 

:U 


2hc^ 


(A.17) 


which  removes  the  spurious  effective  pressure  induced  by  the  FHP  flow’s  kinetic  energy.  The 
ideal  part  of  the  momentum  flux  tensor  takes  the  correct  form 


E 


f 


eq 


bd  2t- 

C  Oij  -\-  pUiUj. 


(A.18) 

(A.19) 


A. 5  Density  Dependent  Pressure 

From  the  PCLBE  in  Section  A. 4  we  know  the  equation  of  state  for  the  isothermal  gas  is  [31] 

P  =  clp.  (A.20) 

We  now  wish  to  consider  how  we  may  alter  the  LBE  to  allow  for  a  more  general  equation  of 
state.  Let  us  add  an  additional  term,  ha,  to  the  R.H.S.  of  (A. 9) 

^  +  ea  ■  V/,  =  -^Ua  -  fa^^^)  +  K{p).  (A.21) 

In  (A.21)  we  have  written  ha  as  a  function  of  the  local  density.  In  a  multi-speed  model 
[33,  34]  ha  may  depend  on  the  local  temperature  as  well. 

We  wish  to  constrain  the  form  of  ha  so  as  not  to  violate  continuity.  We  require 

=  (A.22) 

a 

and  when  V/a  =  0  , 

Eea/ia  =  0.  (A.23) 
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Constraint  (A. 23)  is  required  only  in  uniform  flow;  i.e.  for  general  flows  J2a  ^aha  is  non-zero. 
In  the  uniform  flow  limit  the  LBE  reduces  to 


dfa 

dt 


+  ha{p)- 


(A.24) 


Summing  over  all  lattice  directions  and  using  constraint  (A. 22)  we  have  maintained  the 
collision  property  that 

E^a  =  0.  (A.25) 

a 

Multiplying  by  da,  summing  over  directions,  and  using  (A. 23)  similarly  yields 


^eafla  =  0.  (A.26) 

a 

Thus,  for  arbitrary  flows,  summing  the  LBE  over  all  directions  preserves  continuity 


=  +  (A-27) 

a  a 

=  0,  (A.28) 

where  we  have  used  equations  (A.15,A.16,A.22,  and  A.26). 

In  a  similar  fashion,  we  may  arrive  at  Euler’s  equation 

dtipUi)  +  djiUij)  =  E  Gaihaip),  (A.29) 

a 

where  the  momentum  tensor  is 

^  ^  ^ai^aj  fa'  (A. 30) 

a 

The  R.H.S.  of  (A.29)  imparts  an  effective  density  dependent  pressure.  We  may  expand  ha 
as  follows 

ha  =  h^a^  +  Cajdjh^^^  +  ]^eajeakdjdkh^^'>  +  ^eajCakeaidjdkdih^^'^  H -  (A.31) 

Given  constraint  (A. 22)  and  the  identities  listed  in  the  appendix,  we  immediately  see  that 
and  must  vanish.  So  the  form  of  ha  simplihes  to  a  total  convective-gradient 

ha{p)  =  ea  ■  V  ^yfi'i(p)  +  ^^^ff^eakeaidkdig2{p)^  ,  (A.32) 

where  for  future  convenience  we  have  introduced  gi  =  bh^^^ D  and  g2  =  Shh^^^D^D  +  2). 
The  R.H.S.  of  Euler’s  equation  (A.29)  then  becomes 

E  ^aha{p)  =  V  (^i(p)  +  V^P2(p))  •  (A.33) 


•^tE/a  +  Eea-V/, 

a  a 

— >  dtp  +  V  ■  (pu) 
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Substituting  (A. 18)  and  (A. 33)  into  Euler’s  equation  (A. 29)  gives  the  Navier-Stokes  equation 

dt{pUi)  +  djipUiUj)  =  -di  (c^p  -  gi{p)  -  V‘^g2{p))  ■  (A.34) 

Therefore,  we  have  arrived  at  an  arbitrary  equation  state  dehned  by  functions  gi{p)  and 
5^2  (p)-  The  form  of  the  density  dependent  pressure  follows 

P(P)  =  cV  -  9i{p)  -  V^P2(p)-  (A.35) 


A. 6  Conclusion 

We  have  given  a  new  derivation  of  the  lattice  Boltzmann  equation  starting  from  a  sim¬ 
ple  automaton  Hamiltonian  and  Liouville’s  equations.  We  reviewed  the  LBE  method  and 
illustrated  its  flexibility  in  writing  an  analytical  expression  for  the  system’s  equilibrium  dis¬ 
tribution  to  remove  the  spurious  pressure  kinetic-energy  dependent  term  characteristic  of 
an  FHP  gas.  Now  the  hydrostatic  pressure  correctly  depends  linearly  on  the  local  pressure. 
Given  this  context,  we  generalized  the  LBE  by  the  addition  of  a  convective-gradient  term 
allowing  us  to  model  a  hydrodynamics  governed  by  an  arbitrary  equation  of  state. 


A. 7  Identities  for  Ca 


Our  momentum  states  are  just  the  complex  6*^  roots  of  one 


— 


(A.36) 


where  a  =  1,  2, . . . ,  6.  The  momentum  state-space  has  cardinality  b.  Lattice  summations 
over  odd  powers  of  da  must  vanish  by  symmetry.  The  following  identities,  listed  up  to  the 
fourth  moment,  hold  for  arbitrary  values  of  b  and  spatial  dimension  D  [8] 


=  0 

a 

(A.37) 

^  ^ai^aj  =  71 

a  ^ 

(A.38) 

O 

(A.39) 

(A.40) 
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